/**
 * @file Restrict_Solver.h
 * @brief An Restriction Solver for Multigrid Solver.
 * It can support other multigrid methods.
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-06-09
 * 
 * @copyright Copyright (c) 2021  Linhuo
 */
#ifndef _PRO2_RESTRIC_SOLVER_H_
#define _PRO2_RESTRIC_SOLVER_H_
#include <iostream>
#include <eigen3/Eigen/Core>
#include <cmath>
typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXd;

class Restrict_Solver
{
public:
    virtual VectorXd Restrict(VectorXd _r) = 0;
};

class Restrict_1D_FullWeighting : public Restrict_Solver
{
public:
    VectorXd Restrict(VectorXd _r)
    {
        VectorXd _u_temp;
        int N = (_r.rows() - 1) / 2;
        _u_temp.setZero(N + 1);
        for (int i = 1; i <= N - 1; i++)
        {
            _u_temp(i) = (_r(2 * i - 1) + 2.0 * _r(2 * i) + _r(2 * i + 1)) / 4.0;
        };
        _u_temp(0) = _r(0);
        _u_temp(N) = _r(2 * N);
        return _u_temp;
    };
};

class Restrict_1D_Injection : public Restrict_Solver
{
public:
    VectorXd Restrict(VectorXd _r)
    {
        VectorXd _u_temp;
        int N = (_r.rows() - 1) / 2;
        _u_temp.setZero(N + 1);
        for (int i = 1; i <= N - 1; i++)
        {
            _u_temp(i) = _r(2 * i);
        };
        _u_temp(0) = _r(0);
        _u_temp(N) = _r(2 * N);
        return _u_temp;
    };
};

class Restrict_2D_FullWeighting : public Restrict_Solver
{
public:
    VectorXd Restrict(VectorXd _r)
    {
        VectorXd _u_temp;
        int N = (sqrt(_r.rows()) - 1) / 2;
        _u_temp.setZero((N + 1) * (N + 1));
        double _temp = 0;
        ///ghost cell
        for (int i = 0; i <= N; i++)
        {
            _u_temp(i) = _r(2 * i);
            _u_temp((N + 1) * N + i) = _r((2 * N + 1) * 2 * N + 2 * i);
            _u_temp(i * (N + 1)) = _r(2 * i * (2 * N + 1));
            _u_temp((i + 1) * (N + 1) - 1) = _r((2 * i + 1) * (2 * N + 1) - 1);
        };
        for (int j = 1; j <= N - 1; j++)
        {
            for (int i = 1; i <= N - 1; i++)
            {
                _temp = (1.0 / 16.0) * (_r((2 * j - 1) * (2 * N + 1) + 2 * i - 1) + _r((2 * j - 1) * (2 * N + 1) + 2 * i + 1) + _r((2 * j + 1) * (2 * N + 1) + 2 * i - 1) + _r((2 * j + 1) * (2 * N + 1) + 2 * i + 1));
                _temp += (2.0 / 16.0) * (_r((2 * j - 1) * (2 * N + 1) + 2 * i) + _r((2 * j + 1) * (2 * N + 1) + 2 * i) + _r((2 * j) * (2 * N + 1) + 2 * i - 1) + _r((2 * j) * (2 * N + 1) + 2 * i + 1));
                _temp += (4.0 / 16.0) * (_r((2 * j) * (2 * N + 1) + 2 * i));
                _u_temp(j * (N + 1) + i) = _temp;
            };
        };
        return _u_temp;
    };
};

class Restrict_2D_Injection : public Restrict_Solver
{
public:
    VectorXd Restrict(VectorXd _r)
    {
        VectorXd _u_temp;
        int N = (sqrt(_r.rows()) - 1) / 2;
        _u_temp.setZero((N + 1) * (N + 1));
        double _temp = 0;
        ///ghost cell
        for (int i = 0; i <= N; i++)
        {
            _u_temp(i) = _r(2 * i);
            _u_temp((N + 1) * N + i) = _r((2 * N + 1) * 2 * N + 2 * i);
            _u_temp(i * (N + 1)) = _r(2 * i * (2 * N + 1));
            _u_temp((i + 1) * (N + 1) - 1) = _r((2 * i + 1) * (2 * N + 1) - 1);
        };
        for (int j = 1; j <= N - 1; j++)
        {
            for (int i = 1; i <= N - 1; i++)
            {
                _u_temp(j * (N + 1) + i) = _r((2 * j) * (2 * N + 1) + 2 * i);
            };
        };
        return _u_temp;
    };
};

#else
//DO NOTHING
#endif